# =============================================================================
# April 24 2025
# R code used to produce all results, tables, and figures in the SI Online Appendix J
# Rebecca Cordell
# Unpacking the Role of In-Group Bias in US Public Opinion on Human Rights Violations
# American Journal of Political Science
# https://doi.org/10.7910/DVN/TGAL7M
# =============================================================================

# Clear work environment
rm(list=ls())

# Install Packages
#install.packages("dplyr")
#install.packages("ggplot2")
#install.packages("forcats")
#install.packages("ggpubr")
#install.packages("ggeasy")
#install.packages("lemon")
#install.packages("sandwich")
#install.packages("survey")
#install.packages("lmtest")
#install.packages("remotes")
#remotes::install_version("cregg", version = "0.4.0")

# Required Packages
library("dplyr")
library("ggplot2")
library("forcats")
library("ggpubr")
library("sandwich")
library("survey")
library("lmtest")
library("cregg")
library("ggeasy")
library("lemon")

options(scipen = 999)
options(warn=-1)

# -----------------------------------------------------------------------------
# Read in data
# -----------------------------------------------------------------------------

hr_survey<-read.csv("cordell_ingroupbiashumanrights_data.csv", header=TRUE, stringsAsFactors = FALSE)

# =============================================================================
# Figure J.1: Target (Race) Subgroup Marginal Means
# =============================================================================

# Convert regression variables into factors
factor_vars <- c("perp", "agent", "type", "scope", "targ_nonstate", "targ_race",
                 "targ_relig", "targ_citiz", "frame", "elite", "resp_asian", "resp_black", "resp_hispanic", "resp_middle_east", "resp_white")
hr_survey[factor_vars] <- lapply(hr_survey[factor_vars], as.factor)

# -----------------------------------------------------------------------------
# Calculate subgroup marginal means
# -----------------------------------------------------------------------------

# Model 1
fig_j1_m1_asian <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_asian)
fig_j1_m1_black <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_black)
fig_j1_m1_hispanic <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_hispanic)
fig_j1_m1_middle_east <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_middle_east)
fig_j1_m1_white <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_white)

# Model 2
fig_j1_m2_asian <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_asian)
fig_j1_m2_black <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_black)
fig_j1_m2_hispanic <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_hispanic)
fig_j1_m2_middle_east <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_middle_east)
fig_j1_m2_white <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_white)

# -----------------------------------------------------------------------------
# Prepare Figures
# -----------------------------------------------------------------------------

# Subset targ_race attribute
filter_race <- function(df, race_col) {
  df %>%
    filter(feature == "targ_race", !!sym(race_col) == 1) %>%
    mutate(attrib = "Target (race)")
}
fig_j1_m1_asian <- filter_race(fig_j1_m1_asian, "resp_asian")
fig_j1_m2_asian <- filter_race(fig_j1_m2_asian, "resp_asian")
fig_j1_m1_black <- filter_race(fig_j1_m1_black, "resp_black")
fig_j1_m2_black <- filter_race(fig_j1_m2_black, "resp_black")
fig_j1_m1_hispanic <- filter_race(fig_j1_m1_hispanic, "resp_hispanic")
fig_j1_m2_hispanic <- filter_race(fig_j1_m2_hispanic, "resp_hispanic")
fig_j1_m1_middle_east <- filter_race(fig_j1_m1_middle_east, "resp_middle_east")
fig_j1_m2_middle_east <- filter_race(fig_j1_m2_middle_east, "resp_middle_east")
fig_j1_m1_white <- filter_race(fig_j1_m1_white, "resp_white")
fig_j1_m2_white <- filter_race(fig_j1_m2_white, "resp_white")

# Combine models
fig_j1_asian_all <- rbind(fig_j1_m1_asian, fig_j1_m2_asian)
fig_j1_black_all <- rbind(fig_j1_m1_black, fig_j1_m2_black)
fig_j1_hispanic_all <- rbind(fig_j1_m1_hispanic, fig_j1_m2_hispanic)
fig_j1_middle_east_all <- rbind(fig_j1_m1_middle_east, fig_j1_m2_middle_east)
fig_j1_white_all <- rbind(fig_j1_m1_white, fig_j1_m2_white)

# Create respondents variable (in-groups vs. out-groups)
fig_j1_asian_all <- fig_j1_asian_all %>%
  mutate(Respondents = case_when(
    grepl("asian", level) ~ "In-group",
    level == "Target (race)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j1_black_all <- fig_j1_black_all %>%
  mutate(Respondents = case_when(
    grepl("black", level) ~ "In-group",
    level == "Target (race)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j1_hispanic_all <- fig_j1_hispanic_all %>%
  mutate(Respondents = case_when(
    grepl("hisp", level) ~ "In-group",
    level == "Target (race)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j1_middle_east_all <- fig_j1_middle_east_all %>%
  mutate(Respondents = case_when(
    grepl("middle eastern", level) ~ "In-group",
    level == "Target (race)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j1_white_all <- fig_j1_white_all %>%
  mutate(Respondents = case_when(
    grepl("white", level) ~ "In-group",
    level == "Target (race)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))

# Set factor order
fig_j1_asian_all <- fig_j1_asian_all %>%
  mutate(level = factor(level, levels = c("Target (race)", "asian", "black", "hispanic", "middle eastern", "white")) %>% fct_rev())
fig_j1_black_all <- fig_j1_black_all %>%
  mutate(level = factor(level, levels = c("Target (race)", "black", "asian", "hispanic", "middle eastern", "white")) %>% fct_rev())
fig_j1_hispanic_all <- fig_j1_hispanic_all %>%
  mutate(level = factor(level, levels = c("Target (race)", "hispanic", "asian", "black", "middle eastern", "white")) %>% fct_rev())
fig_j1_middle_east_all <- fig_j1_middle_east_all %>%
  mutate(level = factor(level, levels = c("Target (race)", "middle eastern", "asian", "black", "hispanic", "white")) %>% fct_rev())
fig_j1_white_all <- fig_j1_white_all %>%
  mutate(level = factor(level, levels = c("Target (race)", "white", "asian", "black", "hispanic", "middle eastern")) %>% fct_rev())

# Separate models
fig_j1_m1_asian <- filter(fig_j1_asian_all, outcome == "outcome1")
fig_j1_m2_asian <- filter(fig_j1_asian_all, outcome == "outcome2")
fig_j1_m1_black <- filter(fig_j1_black_all, outcome == "outcome1")
fig_j1_m2_black <- filter(fig_j1_black_all, outcome == "outcome2")
fig_j1_m1_hispanic <- filter(fig_j1_hispanic_all, outcome == "outcome1")
fig_j1_m2_hispanic <- filter(fig_j1_hispanic_all, outcome == "outcome2")
fig_j1_m1_middle_east <- filter(fig_j1_middle_east_all, outcome == "outcome1")
fig_j1_m2_middle_east <- filter(fig_j1_middle_east_all, outcome == "outcome2")
fig_j1_m1_white <- filter(fig_j1_white_all, outcome == "outcome1")
fig_j1_m2_white <- filter(fig_j1_white_all, outcome == "outcome2")

# -----------------------------------------------------------------------------
# Create Figures
# -----------------------------------------------------------------------------

# Create panel (a) Asian respondents - disapproval forced choice
fig_j1_m1_asian_plot <- ggplot(fig_j1_m1_asian, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("White", "Middle eastern", "Hispanic", "Black", "Asian")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.26, 0.74)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(a) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (b) Asian Respondents - disapproval ratings-based
fig_j1_m2_asian_plot <- ggplot(fig_j1_m2_asian, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("White", "Middle eastern", "Hispanic", "Black", "Asian")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.51, 0.89)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(b) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (c) Black respondents - disapproval forced choice
fig_j1_m1_black_plot <- ggplot(fig_j1_m1_black, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("White", "Middle eastern", "Hispanic", "Asian", "Black")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.26, 0.74)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(c) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (d) Black Respondents - disapproval ratings-based
fig_j1_m2_black_plot <- ggplot(fig_j1_m2_black, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("White", "Middle eastern", "Hispanic", "Asian", "Black")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.51, 0.89)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(d) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (e) Hispanic respondents - disapproval forced-choice
fig_j1_m1_hispanic_plot <- ggplot(fig_j1_m1_hispanic, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("White", "Middle eastern", "Black", "Asian", "Hispanic")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.26, 0.74)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(e) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (f) Hispanic respondents - disapproval ratings-based
fig_j1_m2_hispanic_plot <- ggplot(fig_j1_m2_hispanic, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("White", "Middle eastern", "Black", "Asian", "Hispanic")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.51, 0.89)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(f) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (g) Middle Eastern respondents - disapproval forced choice
fig_j1_m1_middle_east_plot <- ggplot(fig_j1_m1_middle_east, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("White", "Hispanic", "Black", "Asian", "Middle eastern")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.26, 0.74)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(g) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (h) Middle Eastern Respondents - disapproval ratings-based
fig_j1_m2_middle_east_plot <- ggplot(fig_j1_m2_middle_east, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("White", "Hispanic", "Black", "Asian", "Middle eastern")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.51, 0.89)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(h) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (i) White respondents - disapproval forced-choice
fig_j1_m1_white_plot <- ggplot(fig_j1_m1_white, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Middle eastern", "Hispanic", "Black", "Asian", "White")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.26, 0.74)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(i) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (j) White respondents - disapproval ratings-based
fig_j1_m2_white_plot <- ggplot(fig_j1_m2_white, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Middle eastern", "Hispanic", "Black", "Asian", "White")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.51, 0.89)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(j) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Print Figure J.1
fig_j1_white_plot<-ggarrange(fig_j1_m1_white_plot, NULL, fig_j1_m2_white_plot,
                             nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j1_white_plot<-annotate_figure(fig_j1_white_plot, top = text_grob("White respondents", size = 14))
fig_j1_black_plot<-ggarrange(fig_j1_m1_black_plot, NULL, fig_j1_m2_black_plot,
                             nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j1_black_plot<-annotate_figure(fig_j1_black_plot, top = text_grob("Black respondents", size = 14))
fig_j1_hispanic_plot<-ggarrange(fig_j1_m1_hispanic_plot, NULL, fig_j1_m2_hispanic_plot,
                            nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j1_hispanic_plot<-annotate_figure(fig_j1_hispanic_plot, top = text_grob("Hispanic respondents", size = 14))
fig_j1_asian_plot<-ggarrange(fig_j1_m1_asian_plot, NULL, fig_j1_m2_asian_plot,
                             nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j1_asian_plot<-annotate_figure(fig_j1_asian_plot, top = text_grob("Asian respondents", size = 14))
fig_j1_middle_east_plot<-ggarrange(fig_j1_m1_middle_east_plot, NULL, fig_j1_m2_middle_east_plot,
                                   nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j1_middle_east_plot<-annotate_figure(fig_j1_middle_east_plot, top = text_grob("Middle eastern respondents", size = 14))
pdf('figure_j1.pdf', width = 8, height = 15)
ggarrange(fig_j1_asian_plot, fig_j1_black_plot, fig_j1_hispanic_plot,  fig_j1_middle_east_plot, fig_j1_white_plot,
          nrow = 5, common.legend = T, legend = "bottom", ncol=1)
dev.off()

# =============================================================================
# Figure J.2: Target (Religion) Subgroup Marginal Means
# =============================================================================

# Create dummy variables
hr_survey$resp_buddhist <- ifelse(hr_survey$resp_relig == "Buddhist", 1, 0)
hr_survey$resp_hindu <- ifelse(hr_survey$resp_relig == "Hindu", 1, 0)
hr_survey$resp_christian <- ifelse(hr_survey$resp_relig == "Christian", 1, 0)
hr_survey$resp_jewish <- ifelse(hr_survey$resp_relig == "Jewish", 1, 0)
hr_survey$resp_muslim <- ifelse(hr_survey$resp_relig == "Muslim", 1, 0)

# Convert regression variables into factors
factor_vars <- c("perp", "agent", "type", "scope", "targ_nonstate", "targ_race",
                 "targ_relig", "targ_citiz", "frame", "elite", "resp_buddhist", "resp_christian", "resp_hindu", "resp_jewish", "resp_muslim")
hr_survey[factor_vars] <- lapply(hr_survey[factor_vars], as.factor)

# -----------------------------------------------------------------------------
# Calculate subgroup marginal means
# -----------------------------------------------------------------------------

# Model 1
fig_j2_m1_buddhist <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_buddhist)
fig_j2_m1_christian <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_christian)
fig_j2_m1_hindu <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_hindu)
fig_j2_m1_jewish <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_jewish)
fig_j2_m1_muslim <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_muslim)

# Model 2
fig_j2_m2_buddhist <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_buddhist)
fig_j2_m2_christian <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_christian)
fig_j2_m2_hindu <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_hindu)
fig_j2_m2_jewish <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_jewish)
fig_j2_m2_muslim <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_muslim)

# -----------------------------------------------------------------------------
# Prepare Figures
# -----------------------------------------------------------------------------

# Subset targ_relig attribute
filter_relig <- function(df, relig_col) {
  df %>%
    filter(feature == "targ_relig", !!sym(relig_col) == 1) %>%
    mutate(attrib = "Target (religion)")
}
fig_j2_m1_buddhist <- filter_relig(fig_j2_m1_buddhist, "resp_buddhist")
fig_j2_m2_buddhist <- filter_relig(fig_j2_m2_buddhist, "resp_buddhist")
fig_j2_m1_christian <- filter_relig(fig_j2_m1_christian, "resp_christian")
fig_j2_m2_christian <- filter_relig(fig_j2_m2_christian, "resp_christian")
fig_j2_m1_hindu <- filter_relig(fig_j2_m1_hindu, "resp_hindu")
fig_j2_m2_hindu <- filter_relig(fig_j2_m2_hindu, "resp_hindu")
fig_j2_m1_jewish <- filter_relig(fig_j2_m1_jewish, "resp_jewish")
fig_j2_m2_jewish <- filter_relig(fig_j2_m2_jewish, "resp_jewish")
fig_j2_m1_muslim <- filter_relig(fig_j2_m1_muslim, "resp_muslim")
fig_j2_m2_muslim <- filter_relig(fig_j2_m2_muslim, "resp_muslim")

# Combine models
fig_j2_buddhist_all <- rbind(fig_j2_m1_buddhist, fig_j2_m2_buddhist)
fig_j2_christian_all <- rbind(fig_j2_m1_christian, fig_j2_m2_christian)
fig_j2_hindu_all <- rbind(fig_j2_m1_hindu, fig_j2_m2_hindu)
fig_j2_jewish_all <- rbind(fig_j2_m1_jewish, fig_j2_m2_jewish)
fig_j2_muslim_all <- rbind(fig_j2_m1_muslim, fig_j2_m2_muslim)

# Create respondents variable (in-groups vs. out-groups)
fig_j2_buddhist_all <- fig_j2_buddhist_all %>%
  mutate(Respondents = case_when(
    grepl("buddhist", level) ~ "In-group",
    level == "Target (religion)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j2_christian_all <- fig_j2_christian_all %>%
  mutate(Respondents = case_when(
    grepl("christian", level) ~ "In-group",
    level == "Target (religion)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j2_hindu_all <- fig_j2_hindu_all %>%
  mutate(Respondents = case_when(
    grepl("hindu", level) ~ "In-group",
    level == "Target (religion)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j2_jewish_all <- fig_j2_jewish_all %>%
  mutate(Respondents = case_when(
    grepl("jewish", level) ~ "In-group",
    level == "Target (religion)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j2_muslim_all <- fig_j2_muslim_all %>%
  mutate(Respondents = case_when(
    grepl("muslim", level) ~ "In-group",
    level == "Target (religion)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))

# Set factor order
fig_j2_buddhist_all <- fig_j2_buddhist_all %>%
  mutate(level = factor(level, levels = c("Target (religion)", "buddhist", "christian", "hindu", "jewish", "muslim")) %>% fct_rev())
fig_j2_christian_all <- fig_j2_christian_all %>%
  mutate(level = factor(level, levels = c("Target (religion)", "christian", "buddhist", "hindu", "jewish", "muslim")) %>% fct_rev())
fig_j2_hindu_all <- fig_j2_hindu_all %>%
  mutate(level = factor(level, levels = c("Target (religion)", "hindu", "buddhist", "christian", "jewish", "muslim")) %>% fct_rev())
fig_j2_jewish_all <- fig_j2_jewish_all %>%
  mutate(level = factor(level, levels = c("Target (religion)", "jewish", "buddhist", "christian", "hindu", "muslim")) %>% fct_rev())
fig_j2_muslim_all <- fig_j2_muslim_all %>%
  mutate(level = factor(level, levels = c("Target (religion)", "muslim", "buddhist", "christian", "hindu", "jewish")) %>% fct_rev())

# Separate models
fig_j2_m1_buddhist <- filter(fig_j2_buddhist_all, outcome == "outcome1")
fig_j2_m2_buddhist <- filter(fig_j2_buddhist_all, outcome == "outcome2")
fig_j2_m1_christian <- filter(fig_j2_christian_all, outcome == "outcome1")
fig_j2_m2_christian <- filter(fig_j2_christian_all, outcome == "outcome2")
fig_j2_m1_hindu <- filter(fig_j2_hindu_all, outcome == "outcome1")
fig_j2_m2_hindu <- filter(fig_j2_hindu_all, outcome == "outcome2")
fig_j2_m1_jewish <- filter(fig_j2_jewish_all, outcome == "outcome1")
fig_j2_m2_jewish <- filter(fig_j2_jewish_all, outcome == "outcome2")
fig_j2_m1_muslim <- filter(fig_j2_muslim_all, outcome == "outcome1")
fig_j2_m2_muslim <- filter(fig_j2_muslim_all, outcome == "outcome2")

# -----------------------------------------------------------------------------
# Create Figures
# -----------------------------------------------------------------------------

# Create panel (a) Buddhist respondents - disapproval forced choice
fig_j2_m1_buddhist_plot <- ggplot(fig_j2_m1_buddhist, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Muslim", "Jewish", "Hindu", "Christian", "Buddhist")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.02, 0.98)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(a) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (b) Buddhist Respondents - disapproval ratings-based
fig_j2_m2_buddhist_plot <- ggplot(fig_j2_m2_buddhist, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Muslim", "Jewish", "Hindu", "Christian", "Buddhist")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.29, 0.96)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(b) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (c) Christian respondents - disapproval forced choice
fig_j2_m1_christian_plot <- ggplot(fig_j2_m1_christian, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Muslim", "Jewish", "Hindu", "Buddhist", "Christian")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.02, 0.98)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(c) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (d) Christian Respondents - disapproval ratings-based
fig_j2_m2_christian_plot <- ggplot(fig_j2_m2_christian, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Muslim", "Jewish", "Hindu", "Buddhist", "Christian")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.29, 0.96)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(d) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (e) Hindu respondents - disapproval forced choice
fig_j2_m1_hindu_plot <- ggplot(fig_j2_m1_hindu, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Muslim", "Jewish", "Christian", "Buddhist", "Hindu")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.02, 0.98)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(e) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (f) Hindu Respondents - disapproval ratings-based
fig_j2_m2_hindu_plot <- ggplot(fig_j2_m2_hindu, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Muslim", "Jewish", "Christian", "Buddhist", "Hindu")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.29, 0.96)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(f) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (g) Jewish respondents - disapproval forced-choice
fig_j2_m1_jewish_plot <- ggplot(fig_j2_m1_jewish, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Muslim", "Hindu", "Christian", "Buddhist", "Jewish")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.02, 0.98)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(g) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (h) Jewish respondents - disapproval ratings-based
fig_j2_m2_jewish_plot <- ggplot(fig_j2_m2_jewish, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Muslim", "Hindu", "Christian", "Buddhist", "Jewish")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.29, 0.96)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(h) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (i) Muslim respondents - disapproval forced-choice
fig_j2_m1_muslim_plot <- ggplot(fig_j2_m1_muslim, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Jewish", "Hindu", "Christian", "Buddhist", "Muslim")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.02, 0.98)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(i) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (j) Muslim respondents - disapproval ratings-based
fig_j2_m2_muslim_plot <- ggplot(fig_j2_m2_muslim, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Jewish", "Hindu", "Christian", "Buddhist", "Muslim")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.29, 0.96)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(j) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Print Figure J.2 
fig_j2_buddhist_plot<-ggarrange(fig_j2_m1_buddhist_plot, NULL, fig_j2_m2_buddhist_plot,
                                nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j2_buddhist_plot<-annotate_figure(fig_j2_buddhist_plot, top = text_grob("Buddhist respondents", size = 14))
fig_j2_christian_plot<-ggarrange(fig_j2_m1_christian_plot, NULL, fig_j2_m2_christian_plot,
                                 nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j2_christian_plot<-annotate_figure(fig_j2_christian_plot, top = text_grob("Christian respondents", size = 14))
fig_j2_hindu_plot<-ggarrange(fig_j2_m1_hindu_plot, NULL, fig_j2_m2_hindu_plot,
                             nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j2_hindu_plot<-annotate_figure(fig_j2_hindu_plot, top = text_grob("Hindu respondents", size = 14))
fig_j2_jewish_plot<-ggarrange(fig_j2_m1_jewish_plot, NULL, fig_j2_m2_jewish_plot,
                              nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j2_jewish_plot<-annotate_figure(fig_j2_jewish_plot, top = text_grob("Jewish respondents", size = 14))
fig_j2_muslim_plot<-ggarrange(fig_j2_m1_muslim_plot, NULL, fig_j2_m2_muslim_plot,
                              nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j2_muslim_plot<-annotate_figure(fig_j2_muslim_plot, top = text_grob("Muslim respondents", size = 14))
pdf('figure_j2.pdf', width = 8, height = 15)
ggarrange(fig_j2_buddhist_plot, fig_j2_christian_plot, fig_j2_hindu_plot, fig_j2_jewish_plot, fig_j2_muslim_plot,
          nrow = 5, common.legend = T, legend = "bottom", ncol=1)
dev.off()

# =============================================================================
# Figure J.3: Target (Citizenship) Subgroup Marginal Means
# =============================================================================

# Create dummy variables
hr_survey$resp_nativ_citiz <- ifelse(hr_survey$resp_citiz %in% c("First generation", "Second generation", "Third generation"), 1, 0)
hr_survey$resp_immig_citiz <- ifelse(hr_survey$resp_citiz == "Immigrant Citizen", 1, 0)
hr_survey$resp_immig_noncitiz <- ifelse(hr_survey$resp_citiz == "Immigrant non-citizen", 1, 0)

# Convert regression variables into factors
factor_vars <- c("perp", "agent", "type", "scope", "targ_nonstate", "targ_race",
                 "targ_relig", "targ_citiz", "frame", "elite", "resp_nativ_citiz", "resp_immig_citiz", "resp_immig_noncitiz")
hr_survey[factor_vars] <- lapply(hr_survey[factor_vars], as.factor)

# -----------------------------------------------------------------------------
# Calculate subgroup marginal means
# -----------------------------------------------------------------------------

# Model 1
fig_j3_m1_nativ_citiz <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_nativ_citiz)
fig_j3_m1_immig_citiz <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_immig_citiz)
fig_j3_m1_immig_noncitiz <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_immig_noncitiz)

# Model 2
fig_j3_m2_nativ_citiz <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_nativ_citiz)
fig_j3_m2_immig_citiz <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_immig_citiz)
fig_j3_m2_immig_noncitiz <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_immig_noncitiz)

# -----------------------------------------------------------------------------
# Prepare Figures
# -----------------------------------------------------------------------------

# Subset targ_citiz attribute
filter_citiz <- function(df, relig_col) {
  df %>%
    filter(feature == "targ_citiz", !!sym(relig_col) == 1) %>%
    mutate(attrib = "Target (citizenship)")
}
fig_j3_m1_nativ_citiz <- filter_citiz(fig_j3_m1_nativ_citiz, "resp_nativ_citiz")
fig_j3_m2_nativ_citiz <- filter_citiz(fig_j3_m2_nativ_citiz, "resp_nativ_citiz")
fig_j3_m1_immig_citiz <- filter_citiz(fig_j3_m1_immig_citiz, "resp_immig_citiz")
fig_j3_m2_immig_citiz <- filter_citiz(fig_j3_m2_immig_citiz, "resp_immig_citiz")
fig_j3_m1_immig_noncitiz <- filter_citiz(fig_j3_m1_immig_noncitiz, "resp_immig_noncitiz")
fig_j3_m2_immig_noncitiz <- filter_citiz(fig_j3_m2_immig_noncitiz, "resp_immig_noncitiz")

# Combine models
fig_j3_nativ_citiz_all <- rbind(fig_j3_m1_nativ_citiz, fig_j3_m2_nativ_citiz)
fig_j3_immig_citiz_all <- rbind(fig_j3_m1_immig_citiz, fig_j3_m2_immig_citiz)
fig_j3_immig_noncitiz_all <- rbind(fig_j3_m1_immig_noncitiz, fig_j3_m2_immig_noncitiz)

# Create respondents variable (in-groups vs. out-groups)
fig_j3_nativ_citiz_all <- fig_j3_nativ_citiz_all %>%
  mutate(Respondents = case_when(
    level == "American citizens" ~ "In-group",
    level == "Target (citizenship)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j3_immig_citiz_all <- fig_j3_immig_citiz_all %>%
  mutate(Respondents = case_when(
    grepl("naturalized American citizens", level) ~ "In-group",
    level == "Target (citizenship)" ~ "In-group",
    TRUE ~ "Out-group"
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j3_immig_noncitiz_all <- fig_j3_immig_noncitiz_all %>%
  mutate(Respondents = case_when(
    level %in% c("immigrants with legal status", "immigrants without legal status") ~ "In-group",
    level == "Target (citizenship)" ~ "In-group",
    TRUE ~ "Out-group"
  ))

# Set factor order
fig_j3_nativ_citiz_all <- fig_j3_nativ_citiz_all %>%
  mutate(level = factor(level, levels = c("Target (citizenship)", "American citizens", "naturalized American citizens", "immigrants with legal status", "immigrants without legal status")) %>% fct_rev())
fig_j3_immig_citiz_all <- fig_j3_immig_citiz_all %>%
  mutate(level = factor(level, levels = c("Target (citizenship)", "naturalized American citizens", "American citizens", "immigrants with legal status", "immigrants without legal status")) %>% fct_rev())
fig_j3_immig_noncitiz_all <- fig_j3_immig_noncitiz_all %>%
  mutate(level = factor(level, levels = c("Target (citizenship)", "immigrants with legal status", "immigrants without legal status", "naturalized American citizens", "American citizens")) %>% fct_rev())

# Separate models
fig_j3_m1_nativ_citiz <- filter(fig_j3_nativ_citiz_all, outcome == "outcome1")
fig_j3_m2_nativ_citiz <- filter(fig_j3_nativ_citiz_all, outcome == "outcome2")
fig_j3_m1_immig_citiz <- filter(fig_j3_immig_citiz_all, outcome == "outcome1")
fig_j3_m2_immig_citiz <- filter(fig_j3_immig_citiz_all, outcome == "outcome2")
fig_j3_m1_immig_noncitiz <- filter(fig_j3_immig_noncitiz_all, outcome == "outcome1")
fig_j3_m2_immig_noncitiz <- filter(fig_j3_immig_noncitiz_all, outcome == "outcome2")

# -----------------------------------------------------------------------------
# Create Figures
# -----------------------------------------------------------------------------

# Create panel (a) American citizens respondents - disapproval forced choice
fig_j3_m1_nativ_citiz_plot <- ggplot(fig_j3_m1_nativ_citiz, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Immigr undoc", "Immigr doc", "Naturalized", "Citizens")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.40, 0.60)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(a) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (b) American citizens Respondents - disapproval ratings-based
fig_j3_m2_nativ_citiz_plot <- ggplot(fig_j3_m2_nativ_citiz, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Immigr undoc", "Immigr doc", "Naturalized", "Citizens")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.48, 0.74)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(b) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (c) naturalized American citizens respondents - disapproval forced choice
fig_j3_m1_immig_citiz_plot <- ggplot(fig_j3_m1_immig_citiz, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Immigr undoc", "Immigr doc", "Citizens", "Naturalized")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.40, 0.60)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(c) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (d) naturalized American citizens Respondents - disapproval ratings-based
fig_j3_m2_immig_citiz_plot <- ggplot(fig_j3_m2_immig_citiz, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Immigr undoc", "Immigr doc", "Citizens", "Naturalized")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.48, 0.74)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(d) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (e) immigrants respondents - disapproval forced choice
fig_j3_m1_immig_noncitiz_plot <- ggplot(fig_j3_m1_immig_noncitiz, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Citizens", "Naturalized", "Immigr undoc", "Immigr doc")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.40, 0.60)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(e) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (f) immigrants Respondents - disapproval ratings-based
fig_j3_m2_immig_noncitiz_plot <- ggplot(fig_j3_m2_immig_noncitiz, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Citizens", "Naturalized", "Immigr undoc", "Immigr doc")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.48, 0.74)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(f) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Print Figure J.3
fig_j3_nativ_citiz_plot<-ggarrange(fig_j3_m1_nativ_citiz_plot, NULL, fig_j3_m2_nativ_citiz_plot,
                                   nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j3_nativ_citiz_plot<-annotate_figure(fig_j3_nativ_citiz_plot, top = text_grob("Native−born citizen respondents", size = 14))
fig_j3_immig_citiz_plot<-ggarrange(fig_j3_m1_immig_citiz_plot, NULL, fig_j3_m2_immig_citiz_plot,
                                   nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j3_immig_citiz_plot<-annotate_figure(fig_j3_immig_citiz_plot, top = text_grob("Immigrant citizen respondents", size = 14))
fig_j3_immig_noncitiz_plot<-ggarrange(fig_j3_m1_immig_noncitiz_plot, NULL, fig_j3_m2_immig_noncitiz_plot,
                                      nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j3_immig_noncitiz_plot<-annotate_figure(fig_j3_immig_noncitiz_plot, top = text_grob("Immigrant non−citizen respondents", size = 14))
pdf('figure_j3.pdf', width = 8, height = 11)
ggarrange(fig_j3_nativ_citiz_plot, fig_j3_immig_citiz_plot, fig_j3_immig_noncitiz_plot,
          nrow = 3, common.legend = T, legend = "bottom", ncol=1)
dev.off()

# =============================================================================
# Figure J.4: Elite Cue (Partisanship) Subgroup Marginal Means
# =============================================================================

# Convert regression variables into factors
factor_vars <- c("perp", "agent", "type", "scope", "targ_nonstate", "targ_race",
                 "targ_relig", "targ_citiz", "frame", "elite", "resp_dem", "resp_rep")
hr_survey[factor_vars] <- lapply(hr_survey[factor_vars], as.factor)

# Model 1
fig_j4_m1_dem <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_dem)
fig_j4_m1_rep <- cregg::cj(hr_survey, outcome1 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_rep)

# Model 2
fig_j4_m2_dem <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_dem)
fig_j4_m2_rep <- cregg::cj(hr_survey, outcome2 ~ perp + agent + type + scope + targ_nonstate + targ_race + targ_relig + targ_citiz + frame + elite, id = ~id, estimate = "mm", by = ~resp_rep)

# -----------------------------------------------------------------------------
# Prepare Figures
# -----------------------------------------------------------------------------

# Subset elite attribute
filter_elite <- function(df, party_col) {
  df %>%
    filter(
      grepl("A Democrat member of Congress|A Republican member of Congress", level),
      feature == "elite",
      !!sym(party_col) == 1
    ) %>%
    mutate(attrib = "Elite (Partisanship)")
}
fig_j4_m1_dem <- filter_elite(fig_j4_m1_dem, "resp_dem")
fig_j4_m2_dem <- filter_elite(fig_j4_m2_dem, "resp_dem")
fig_j4_m1_rep <- filter_elite(fig_j4_m1_rep, "resp_rep")
fig_j4_m2_rep <- filter_elite(fig_j4_m2_rep, "resp_rep")

# Combine models
fig_j4_dem_all <- rbind(fig_j4_m1_dem, fig_j4_m2_dem)
fig_j4_rep_all <- rbind(fig_j4_m1_rep, fig_j4_m2_rep)

# Create respondents variable (in-groups vs. out-groups)
fig_j4_dem_all <- fig_j4_dem_all %>%
  mutate(Respondents = case_when(
    grepl("Democrat", level) ~ "In-group",
    grepl("Republican", level) ~ "Out-group",
    level == "Elite (Partisanship)" ~ "In-group",
    TRUE ~ level
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))
fig_j4_rep_all <- fig_j4_rep_all %>%
  mutate(Respondents = case_when(
    grepl("Republican", level) ~ "In-group",
    grepl("Democrat", level) ~ "Out-group",
    level == "Elite (Partisanship)" ~ "In-group",
    TRUE ~ level
  )) %>%
  mutate(Respondents = factor(Respondents, levels = c("In-group", "Out-group")))

# Set factor order
fig_j4_dem_all <- fig_j4_dem_all %>%
  mutate(level = factor(level, levels = c("Elite (Partisanship)", "A Democrat member of Congress", "A Republican member of Congress")) %>% fct_rev())
fig_j4_rep_all <- fig_j4_rep_all %>%
  mutate(level = factor(level, levels = c("Elite (Partisanship)", "A Republican member of Congress", "A Democrat member of Congress")) %>% fct_rev())

# Separate models
fig_j4_m1_dem <- filter(fig_j4_dem_all, outcome == "outcome1")
fig_j4_m2_dem <- filter(fig_j4_dem_all, outcome == "outcome2")
fig_j4_m1_rep <- filter(fig_j4_rep_all, outcome == "outcome1")
fig_j4_m2_rep <- filter(fig_j4_rep_all, outcome == "outcome2")

# -----------------------------------------------------------------------------
# Create Figures
# -----------------------------------------------------------------------------

# Create panel (a) Democrat respondents - disapproval forced-choice
fig_j4_m1_dem_plot <- ggplot(fig_j4_m1_dem, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Republican MOC", "Democrat MOC")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.46, 0.54)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(a) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (b) Democrat respondents - disapproval ratings-based
fig_j4_m2_dem_plot <- ggplot(fig_j4_m2_dem, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Republican MOC", "Democrat MOC")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.62, 0.76)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(b) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  labs(color = NULL)

# Create panel (c) Republican respondents - disapproval forced choice
fig_j4_m1_rep_plot <- ggplot(fig_j4_m1_rep, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Democrat MOC", "Republican MOC")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.46, 0.54)) +
  geom_hline(yintercept = 0.5, size = 0.2, color = "black") +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(c) Disapproval forced-choice") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Create panel (d) Republican Respondents - disapproval ratings-based
fig_j4_m2_rep_plot <- ggplot(fig_j4_m2_rep, aes(x = level, y = estimate, colour = Respondents)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), size = 0.3, width = 0.2) +
  scale_x_discrete(labels = c("Democrat MOC", "Republican MOC")) +
  xlab('') + ylab('Marginal mean') +
  coord_flip(ylim = c(0.62, 0.76)) +
  theme_classic(base_size = 10) +
  scale_colour_grey() +
  easy_center_title() +
  ggtitle("(d) Disapproval ratings-based") +
  theme(
    axis.text.x = element_text(colour = "black"),
    axis.text.y = element_text(colour = "black"),
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.ticks.y = element_blank()
  ) +
  labs(color = NULL)

# Print Figure J.4
fig_j4_dem_plot<-ggarrange(fig_j4_m1_dem_plot, NULL, fig_j4_m2_dem_plot,
                           nrow = 1, common.legend = F, legend = NULL, ncol = 3, widths = c(1, 0.09, 1))
fig_j4_dem_plot<-annotate_figure(fig_j4_dem_plot, top = text_grob("Democrat respondents", size = 14))
fig_j4_rep_plot<-ggarrange(fig_j4_m1_rep_plot, NULL, fig_j4_m2_rep_plot,
                           nrow = 1, common.legend = T, legend = "bottom", ncol = 3, widths = c(1, 0.09, 1))
fig_j4_rep_plot<-annotate_figure(fig_j4_rep_plot, top = text_grob("Republican respondents", size = 14))
pdf('figure_j4.pdf', width = 8.26, height = 5.82)
ggarrange(fig_j4_dem_plot, fig_j4_rep_plot,
          nrow = 2, common.legend = T, legend = "bottom", ncol=1, heights = c(1, 1.1))
dev.off()